#include "Random.hpp"
#include <armadillo>
#include <boost/math/distributions/normal.hpp>
#include "utils.hpp"
#include <boost/math/distributions/students_t.hpp>
#include <boost/random/student_t_distribution.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <cmath> 
#include "simulate/Owen.hpp"
#include <boost/math/distributions/gamma.hpp>
#include <boost/random/binomial_distribution.hpp>
//#include <boost/math/distribution/binomial.hpp>
#include <boost/math/distributions/binomial.hpp>

namespace Distribution{
#define _mINF -100000000000
#define tZERO 0.000000000001

using namespace arma;

class Distribution
{
	public:
	virtual mat r(int)=0;
	virtual double d(mat,int,int)=0;
	virtual mat Var(mat)=0;
	virtual mat scrambled(int)=0;
	virtual mat Get_s()=0;
	virtual mat Get_m()=0;
};

class Gaussian : public Distribution
{
	public:
	Gaussian(){}
	Gaussian(int n,mat m, mat s)
	{
		_n=n;
		_m=m;
		try{
			_s=s;//chol here
			_cs=chol(s);
			_is=inv(s);
		}catch(std::exception e)
		{
			_s=diagmat(s);
			_cs=sqrt(_s);
			_is=1/(_s);
			cout << "big problem";
			
		}
		boost::random::normal_distribution<> N;
		dn= new boost::math::normal_distribution<>(0,1); 
		rnorm= new RandomG::Random<boost::random::normal_distribution<> >(N);
	}
	virtual ~Gaussian()
	{
		delete rnorm;
		delete dn;
	}

	mat r(int j)
	{
		mat X(_n,1);
		for(int i=0;i<_n;i++)
		{
		 	X(i,0)=(*rnorm)();
		}
	 	X=_m+(_cs).t()*X;
		return X;
	}


	double 	d(mat x,int blog=0,int kk=0){
	
if(blog){
		double res=0;
		mat foo=x-_m;
		double bar=as_scalar(foo.t()*_is*foo);
		//cout << "bar" << bar;
		res=-_n*0.5*log(2*PI)-0.5*log(det(_s))-0.5*bar;
	//	cout << "res" << log(det(_s));
		return res;

}else{
		double res=0;
		mat foo=x-_m;
		mat bar=foo.t()*_is*foo;
		double bar1=as_scalar(bar);
		long double temp=det(_s);
		long double st=sqrt(temp);
		long double r=(sqrt(pow(2*PI,_n))*st);
		res=(double)1/r;
		double res2=res*exp(-0.5*bar1);
		if((res2>_mINF) & (res2!=0))
		{
			return res2;
		}else if(res2==0){
			return tZERO;
		}else if(res2<_mINF){
			return _mINF;
		}
}
	}
	mat scrambled(int M)
	{
		int Mn=M*_n+1;
		double *N=new double[Mn];
		sobol::Owenscramble O;
		O(M,_n,1,N);
		mat X=vect2mat<double>(N,M,_n);
		delete[] N;
		mat temp;
		for( int i=0;i<M;i++)
		{
			temp=(_m+_cs.t()*X(i,span::all).t());
		//	temp=this->r(1);
			X(i,span::all)=temp.t();
		}
		return X;
	}
	mat Var(mat X){return arma::cov(X);}
	void Set_n(int n){	_n=n;	}
	void Set_m(mat m){	_m=m;	}
	void Set_s(mat s){
		_s=s;
		try{
			_cs=chol(s);
			_is=inv(s);
		}catch(std::exception e)
		{
			cout << "Variance covariance is singular";
			_s=diagmat(s);
			_cs=sqrt(_s);
			_is=inv(_s);
		}
	}
	mat Get_s(void){return _s;}
	mat Get_m(void){return _m;}
	private:
	int _n;
	mat _m;
	mat _s;
	mat _is;
	mat _cs;
	RandomG::Random<boost::random::normal_distribution<> > *rnorm;
	boost::math::normal_distribution<> *dn;

};

#include "Student.hpp"
#include "Uniform.hpp"
#include "LogitProp.hpp"
#include "InvKhi2.hpp"
#include "Binomial.hpp"
}
